function G = slnngraph(X, X2, k, varargin)
%SLNNGRAPH Creates a nearest neighbor graph
%
% [ Syntax ]
%   - G = slnngraph(X, [], k, ...)
%   - G = slnngraph(X, X, k, ...)
%   - G = slnngraph(X, X2, k, ...)
%
% [ Arguments ]
%   - X:        the source sample set
%   - X2:       the destination sample set
%   - k:        the (maximum) number of neighbors for each sample
%   - G:        the created graph (in adjacent matrix)
%
% [ Description ]
%   - G = slnngraph(X, [], k, ...) creates a nearest neighbor graph by
%     connecting each node and its neighbors. 
%
%     Suppose there are n samples in d-dimensional space, then X should be
%     a d x n matrix. The created graph has n nodes, each node links to its
%     k nearest neighbors (the node itself is excluded from its neighbor
%     set). Correspondingly, G is a n x n full or sparse matrix (by
%     default, it is sparse).
%
%     If j-th node is one of the neighbors of i-th node, then G(i, j) is
%     the value associated with the edge linking i-th node and j-th
%     node; otherwise, G(i, j) is zero. By default, the value for each edge
%     is set to 1. However, we can compute the edge values based on the
%     distance between the nodes that is links or the nodes themselves by
%     specifying the dfun or efun options.
%
%   - G = slnngraph(X, X, k, ...) creates a nearest neighbor graph by
%     connecting each node and its neighbors. By this syntax, each node
%     itself is counted as one of its own neighbors.
%
%   - G = slnngraph(X, X2, k, ...) creates a nearest neighbor bigraph by
%     connecting each node in X to its nearest neighbors in X2. 
%
%     Here, X and X2 can have the same number or different numbers of 
%     samples. Suppose X and X2 respectively have n and n2 samples. Then
%     the created graph G will be an n x n2 matrix. 
%
%     You can specify the following options for graph construction in form 
%     of name-value pairs.
%     \{:
%         - output:         the form of output adjacency matrix
%                           (default = 'sparse')
%                           \{:
%                               - full:     G will be a full matrix
%                               - sparse:   G will be a sparse matrix
%                           \:}
%
%         - method:         The method to find nearest neighbors. 
%                           (default = 'std')
%                           \{:
%                               - std:      use function slfindnn, which 
%                                           searches nearest neighbors by
%                                           computing all pairwise
%                                           distances.
%                               - ann:      use annquery in the ANN MATLAB
%                                           Wrapper, which is another
%                                           toolbox available in MATLAB
%                                           Central File Exchange.
%
%                                           This method can only be used
%                                           when ANN MATLAB Wrapper is
%                                           installed.
%                           \:}
%
%         - metric:         The distance metric, which can be either a
%                           metric name available in function slmetric_pw,
%                           or a pairwise metric function handle that
%                           supports the following syntax
%                              $ D = f(X1, X2) $
%                           where X1 and X2 are d x n1 and d x n2 sample
%                           matrices. Then D should be a n1 x n2 matrix
%                           of pairwise metric values.
%                           By default, metric is 'eucdist'.
%
%                           This option only take effects for std method.
%                           The ann method can only use euclidean distance.
%
%         - dfun:           The function handle to compute edge values based 
%                           on the distance metrics between two end nodes.
%                           It should support the following syntax
%                               $ v = f(mv) $
%                           Here, mv is a column vector of metric values, v 
%                           should be a vector of the same size, giving the
%                           edge values.
%
%                           By default, it is [], which means dfun is not
%                           specified.
%
%         - efun:           The function handle to compute edge values based 
%                           on the nodes that link. It should support the
%                           following syntax
%                               $ v = f(Y1, Y2) $
%                           Here, Y1 and Y2 are d x ns matrices, which
%                           respectively contain ns samples. v should be 
%                           a 1 x ns row vector, with v(i) corresponding to
%                           the edge linking the samples Y1(:, i) and Y2(:, i).
%                           
%                           By default, it is [], which means efun is not
%                           specified.
%
%         - eval:           The default edge value. (default = 1).
%                           This is the values set to each edge when both
%                           dfun and efun are not specified.
%
%         - radius:         The upper bound on the distances between neighborinig 
%                           samples.  
%
%                           By default, radius is 0, which means no constraint is 
%                           set to the neighboring distances.
%
%                           If radius is set to a positive value, then only the 
%                           nodes whose distance is not larger than radius
%                           can be connected with an edge.
%
%         - dcache:         The maximum number of distance values in cache.
%                           (default = 1e7). This option only takes effect
%                           when using method std.
%
%         - annoptions:     The cell array of name value pairs to specify
%                           the options for ANN search. This option only
%                           takes effect when using method ann.
%                           default = {}.
%     \:}
%
% [ Remarks ]
%   - The edge values are determined according to the following rules
%       - If efun is specified, then use efun to compute edge values based
%         on the nodes that the edges link.
%       - If efun is not specified and dfun is specified, then use dfun to
%         compute the edge values based on the corresponding distance
%         metric values.
%       - If both efun and dfun is not specified, then the value for each
%         edge is set to eval, whose default value is 1.
%
% [ History ]
%   - Created by Dahua Lin, on Jul 7, 2007
%

%% parse and verify input arguments

error(nargchk(3, inf, nargin));

% X and X2
assert(isnumeric(X) && ndims(X) == 2, 'sltoolbox:slnngraph:invalidarg', ...
    'The X should be a 2D numeric matrix.');

if isempty(X2)
    X2 = X;
    exclude_self = true;
else
    assert(isnumeric(X2) && ndims(X2) == 2, 'sltoolbox:slnngraph:invalidarg', ...
        'The X2 should be a 2D numeric matrix.');
    assert(size(X, 1) == size(X2, 1), 'sltoolbox:slnngraph:dimmismatch', ...
        'X and X2 should be with the same sample dimension.');        
    exclude_self = false;
end
n = size(X, 2);
n2 = size(X2, 2);


% k
assert(isnumeric(k) && isscalar(k) && k > 0 && k == fix(k), 'sltoolbox:slnngraph:invalidarg', ...
    'k should be a positive integer scalar.');

if exclude_self
    assert(k < n, 'sltoolbox:slnngraph:toolarge_k', ...
        'The value of k should be less than the number of nodes in self-exclusion mode.');
else
    assert(k <= n2, 'sltoolbox:slnngraph:toolarge_k', ...
        'The value of k should not be larger than the number of target nodes.');
end

% options

opts = struct( ...
    'output', 'sparse', ...
    'method', 'std', ...
    'metric', 'eucdist', ...
    'dfun', [], ...
    'efun', [], ...
    'eval', 1, ...
    'radius', 0, ...
    'dcache', 1e7, ...
    'annoptions', {{}});

if ~isempty(varargin)
    opts = setopts(opts, varargin{:});
end

% output
assert(ischar(opts.output), 'sltoolbox:slnngraph:invalidopt', ...
    'The option output should be a string');
out_sparse = find(strcmp(opts.output, {'sparse', 'full'})) == 1;
assert(~isempty(out_sparse), 'sltoolbox:slnngraph:invalidopt', ...
    'The option output should be either ''sparse'' or ''full''');

% method
assert(ischar(opts.method), 'sltoolbox:slnngraph:invalidopt', ...
    'The option method should be a string');
use_ann = find(strcmp(opts.method, {'std', 'ann'})) == 2;
assert(~isempty(use_ann), 'sltoolbox:slnngraph:invalidopt', ...
    'The option method should be either ''std'' or ''ann''');
if use_ann
    assert(exist('annquery', 'file') == 2, 'sltoolbox:slnngraph:ann_notfound', ...
        'The ANN MATLAB Wrapper is not installed. The method ann cannot be used.');
end

% metric
if ~use_ann
    assert(ischar(opts.metric) || isa(opts.metric, 'function_handle'), ...
        'sltoolbox:slnngraph:invalidopt', 'The option metric should be either a name or a function handle.');
end

% dun
assert(isempty(opts.dfun) || isa(opts.dfun, 'function_handle'), ...
    'sltoolbox:slnngraph:invalidopt', 'The option dfun should be either empty or a function handle.');
% efun
assert(isempty(opts.efun) || isa(opts.efun, 'function_handle'), ...
    'sltoolbox:slnngraph:invalidopt', 'The option efun should be either empty or a function handle.');
% eval
assert(isnumeric(opts.eval) && isscalar(opts.eval), ...
    'sltoolbox:slnngraph:invalidopt', 'The option eval should be a numeric scalar.');
% radius
assert(isnumeric(opts.radius) && isscalar(opts.radius) && opts.radius >= 0, ...
    'sltoolbox:slnngraph:invalidopt', 'The option radius should be a nonnegative numeric scalar.');
% dcache
if ~use_ann
    assert(isnumeric(opts.dcache) && isscalar(opts.dcache) && opts.dcache > 0, ...
        'sltoolbox:slnngraph:invalidopt', 'The option dcache should be a positive numeric scalar.');
end
% annoptions
if use_ann
    assert(iscell(opts.annoptions), 'sltoolbox:slnngraph:invalidopt', ...
        'The annoptions should be a cell array.');
end

%% main

% find nn and edges

if isempty(opts.efun) && ~isempty(opts.dfun);
    [edges, dists] = find_nn_edges(X, X2, k, exclude_self, use_ann, opts);
else
    edges = find_nn_edges(X, X2, k, exclude_self, use_ann, opts);
end

% allocate graph
if out_sparse
    G = spalloc(n, n2, size(edges, 1));
else
    G = zeros(n, n2);
end

% set edge values
einds = sub2ind([n, n2], edges(:,1), edges(:,2));
if ~isempty(opts.efun)
    evals = opts.efun(X(:, edges(:,1)), X2(:, edges(:,2)));
    G(einds) = evals(:);
    
elseif ~isempty(opts.dfun)
    evals = opts.dfun(dists);
    G(einds) = evals;
    
else
    G(einds) = opts.eval;
end


%% The sub-function to find edges

function [edges, dists] = find_nn_edges(X, X2, k, exclude_self, use_ann, opts)

if ~use_ann
    if ~exclude_self
        Xr = X2;
        Xq = X;
    else
        Xr = X;
        Xq = [];
    end        
    oplist = {'metric', opts.metric, 'radius', opts.radius, 'dcache', opts.dcache};
    
    if nargout < 2
        nnidx = slfindnn(Xr, Xq, k, oplist{:});
        edges = slnn2edges(nnidx);
    else
        [nnidx, dists] = slfindnn(Xr, Xq, k, oplist{:});
        [edges, dists] = slnn2edges(nnidx, dists);
    end
        
else
    Xr = X2;
    Xq = X;
    
    if opts.radius > 0
        oplist = [{'search_sch', 'fr', 'eps', opts.radius}, opts.annoptions];
    else
        oplist = opts.annoptions;
    end
    
    if nargout < 2
        nnidx = annquery(Xr, Xq, k, oplist{:});
        edges = slnn2edges(nnidx);
    else
        [nnidx, dists] = annquery(Xr, Xq, k, oplist{:});
        [edges, dists] = slnn2edges(nnidx, dists);
    end
    
end

 
